//system
#include <sstream>
#include <fstream>
#include <vector>
#include <stdexcept>
#include <deque>
#include <string>
#include <iostream>


//boost
#include <boost/tuple/tuple.hpp>

// my headers
#define DEBUG_LEVEL 6
#include <debug_func.h>
#include <blast_results.h>
#include <index_string.h>
#include <portable_timer.h>
#include <codeml_results.h>
#include <progress_indicator.h>
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <ultoa.h>
#include <print_error.h>
#include <open_fstream.h>
#include <seq_pairs.h>
#include <get_filename.h>
#include <bcpp_read_seq.h>



#include "get_arg.h"


using std::vector;
using std::pair;
using std::deque;
using std::ofstream;
using std::ostream;
using std::ifstream;
using std::istream;
using std::istringstream;
using std::string;
using std::runtime_error;
using std::cerr;
using std::cin;
using std::cout;
using boost::tie;
using BioCPP::t_collect_sequences;
using BioCPP::t_collect_sequences_report_errors_to_strm;
using BioCPP::eSEQ_NODUPLICATES;
using BioCPP::eSEQ_DUPLICATES;
using bcpp_align::t_align_results;



struct t_threshold_counts
{
	string					name;
	unsigned				input;
	unsigned            	ds_threshold;
	unsigned            	dn_threshold;
	unsigned            	tau_threshold;
	unsigned            	len_threshold;
	unsigned            	fragment_coverage_threshold;
	unsigned            	residues_coverage_threshold;
	unsigned            	coverage_threshold;
	unsigned            	codeml_failed;
	unsigned            	output;
	t_program_args&			args;
	t_threshold_counts(const string& name_, t_program_args& args_):
							name(name_),
							input(0),
							ds_threshold(0),
							dn_threshold(0),
							tau_threshold(0),
							len_threshold(0),
							fragment_coverage_threshold(0),
							residues_coverage_threshold(0),
							coverage_threshold(0),
							codeml_failed(0),
							output(0),
							args(args_)
							{}
};

enum E_THRESHOLD_GENE_COUNTS{	eINPUT = 1,
								eDS = 2,
								eDN = 4,
								eTAU = 8,
								eLEN = 16,
								eFRAG = 32,
								eRES = 64,
								eCOVER = 128,
								eCODEML_FAILED = 256,
								eOUTPUT = 512};


string left_pad(unsigned num, unsigned len = 8)
{
	string res = ultoa(num);
	if (res.length() < len)
		res = string(len - res.length(), ' ')  + res;
	return res;
}
ostream& operator << (ostream& os, const t_threshold_counts& counts)
{
	os << "\t"		<< left_pad(counts.input)
								<< " " << counts.name <<  " parsed.\n";
	os << "\t"	<< left_pad(counts.codeml_failed)
							<< " with no codeml output (dNdS estimate failed).\n ";

	if (counts.len_threshold != std::numeric_limits<unsigned>::max())
		os << "\t"	<< left_pad(counts.len_threshold)
								<< " with aligned length >= "
								<< counts.args.len_threshold << "\n";
	if (counts.fragment_coverage_threshold > 0.0)
		os << "\t" 	<< left_pad(counts.fragment_coverage_threshold)
								<< " with fraction of residues aligned (incl. fragments) > "
								<< counts.args.fragment_coverage << "\n";
	if (counts.residues_coverage_threshold > 0.0)
		os << "\t" 	<< left_pad(counts.residues_coverage_threshold)
								<< " with fraction of residues aligned (sans fragments) > "
								<< counts.args.residues_coverage << "\n";
	if (counts.fragment_coverage_threshold > 0.0 || counts.residues_coverage_threshold > 0.0)
		os << "\t" 	<< left_pad(counts.coverage_threshold)
								<< " with coverage above threshold\n";
	if (counts.ds_threshold != std::numeric_limits<float>::max())
		os << "\t" 	<< left_pad(counts.ds_threshold)
								<< " with dS < "
								<< counts.args.ds_threshold << "\n";
	if (counts.dn_threshold != std::numeric_limits<float>::max())
		os << "\t" 	<< left_pad(counts.dn_threshold)
								<< " with dN < "
								<< counts.args.dn_threshold << "\n";
	if (counts.tau_threshold != std::numeric_limits<float>::max())
		os << "\t" 	<< left_pad(counts.tau_threshold)
								<< " with tau < "
								<< counts.args.tau_threshold << "\n";

	os	 << "\t" 	<< left_pad(counts.output)
								<< " " << counts.name <<  " output.\n";

	return os;
}

//________________________________________________________________________________________

//		count_pairs_over_threshold

//________________________________________________________________________________________
void count_pairs_over_threshold(t_threshold_counts&				cnts,
								const std::deque<unsigned>&		gene_cnts)
{
	for (unsigned i = 0; i < gene_cnts.size(); ++i)
	{
		cnts.input							+= gene_cnts[i] & eINPUT;
		cnts.len_threshold					+= gene_cnts[i] & eLEN;
		cnts.ds_threshold					+= gene_cnts[i] & eDS;
		cnts.dn_threshold					+= gene_cnts[i] & eDN;
		cnts.tau_threshold					+= gene_cnts[i] & eTAU;
		cnts.fragment_coverage_threshold	+= gene_cnts[i] & eFRAG;
		cnts.residues_coverage_threshold	+= gene_cnts[i] & eRES;
		cnts.coverage_threshold				+= gene_cnts[i] & eCOVER;
		cnts.output							+= gene_cnts[i] & eOUTPUT;
		cnts.codeml_failed					+= gene_cnts[i] & eCODEML_FAILED;
	}
	cnts.input							/= eINPUT;
	cnts.len_threshold					/= eLEN;
	cnts.ds_threshold					/= eDS;
	cnts.dn_threshold					/= eDN;
	cnts.tau_threshold					/= eTAU;
	cnts.fragment_coverage_threshold	/= eFRAG;
	cnts.residues_coverage_threshold	/= eRES;
	cnts.coverage_threshold				/= eCOVER;
	cnts.output							/= eOUTPUT;
	cnts.codeml_failed					/= eCODEML_FAILED;
}



//________________________________________________________________________________________

//		index_gene_pair

//________________________________________________________________________________________
pair<unsigned, unsigned> index_gene_pair(const string&		name1, 
										 const string&		name2, 
										 t_index_string& 	genes,
										 deque<unsigned>&	gene_cnts)
{
	bool added;
	unsigned gene_index1 = genes.index (name1, added);
	if (added)
		gene_cnts.push_back(1);
	else
		gene_cnts[gene_index1] |= eINPUT;
	unsigned gene_index2 = genes.index (name2, added);
	if (added)
		gene_cnts.push_back(1);
	else
		gene_cnts[gene_index2] |= eINPUT;
	return std::make_pair(gene_index1, gene_index2);
}


//________________________________________________________________________________________

//		do_parse_pairs

//________________________________________________________________________________________
struct do_parse_pairs
{
	private:
		string			old_name1;
		string			old_name2;
		t_index_string	genes;
		const unsigned	coverage_required;
		enum eOUTPUT_MODE {edS, edN, edS_COVERAGE, eALL} output_mode;
	public:
		deque<unsigned>					gene_flags;
		t_threshold_counts&				pair_cnts;
		t_program_args&					args;
		std::string						comment_letters;
		do_parse_pairs(	t_threshold_counts&				pair_cnts_,
						t_program_args&					args_,
						const std::string&				comment_letters_):

				coverage_required(args_.coverage_both ?
												2 :					// both in pair
												1 ),				// either one of pair
				pair_cnts(pair_cnts_),
				args(args_),
				comment_letters(comment_letters_)
		{
			if (args.output_ds_only)
				output_mode = edS;
			else if (args.output_dn_only)
				output_mode = edN;
			else if (args.output_ds_coverage)
				output_mode = edS_COVERAGE;
			else
				output_mode = eALL;
		}
	void print_error(unsigned cnt_line, const string& prot_name) const
	{
		std_print_error(args.CERR(), "The sequence for \"" + prot_name +
								"\" on line #" + ultoa(cnt_line) +
								" could not be found");
	}
	void operator()(istream& is, unsigned cnt_line)
	{
		// get data for one protein pair
		string name1, name2;

		// comment if false
		string comment_line;
		if (!seq_pairs::parse_ids(is, name1, name2, comment_line, comment_letters,
								  seq_pairs::ePARSE_IDS_ONLY))
		{						  
			cout	<< comment_line;
			return;
		}		


		t_codeml_results dndsdata;
		is >> dndsdata;

		t_align_results aligndata;
		is >> aligndata;


		if (!is)
			throw std::runtime_error("The alignment data is not in the correct format");

		
		// prevent simple repeats
		if (name1 == old_name1 && name2 == old_name2)
			return;
		old_name1 = name1;
		old_name2 = name2;

		// index genes
		unsigned gi1, gi2;
		tie(gi1, gi2) = index_gene_pair(name1, name2, genes, gene_flags);
		++pair_cnts.input;

		// check for codeml failures
		if (dndsdata.S == 0.0)
		{
			gene_flags[gi1] |= eCODEML_FAILED;
			gene_flags[gi2] |= eCODEML_FAILED;
			++pair_cnts.codeml_failed;
			return;
		}


		// len threshold
		if (aligndata.align_len < args.len_threshold)
			return;
		{
			gene_flags[gi1] |= eLEN;
			gene_flags[gi2] |= eLEN;
			++pair_cnts.len_threshold;
		}



		//
		// 	Coverage threshold
		// 
		unsigned index_of_fragment = 0;
		{
			unsigned pass_threshold = 0;
			
			// which is shorter
			double min_len, max_len;
			if (aligndata.len1 <= aligndata.len2)
			{
				min_len = aligndata.len1;
				max_len = aligndata.len2;
				index_of_fragment = 1;
			}
			else
			{
				min_len = aligndata.len2;
				max_len = aligndata.len1;
				index_of_fragment = 2;
			}

			
			// fragments threshold
			if (aligndata.aligned_residues >= args.fragment_coverage * min_len)
			{
				gene_flags[gi1] |= eFRAG;
				gene_flags[gi2] |= eFRAG;
				++pair_cnts.fragment_coverage_threshold;
				++pass_threshold;
			}

			// residues threshold
			if (aligndata.aligned_residues >= args.residues_coverage * max_len)
			{
				gene_flags[gi1] |= eRES;
				gene_flags[gi2] |= eRES;
				++pair_cnts.residues_coverage_threshold;
				++pass_threshold;
				// these are not fragments. Set flag to 0
				index_of_fragment = 0;
			}

			// check if pass coverage threshold (and / or)
			if (pass_threshold < coverage_required )
				return;

			gene_flags[gi1] |= eCOVER;
			gene_flags[gi2] |= eCOVER;
			++pair_cnts.coverage_threshold;
		}


		// dS threshold
		if (dndsdata.dS > args.ds_threshold)
			return;
		{
			gene_flags[gi1] |= eDS;
			gene_flags[gi2] |= eDS;
			++pair_cnts.ds_threshold;
		}


		// dN threshold
		if (dndsdata.dN > args.dn_threshold)
			return;
		{
			gene_flags[gi1] |= eDN;
			gene_flags[gi2] |= eDN;
			++pair_cnts.dn_threshold;
		}

		// tau threshold
		if (dndsdata.dT > args.tau_threshold)
			return;
		{
			gene_flags[gi1] |= eTAU;
			gene_flags[gi2] |= eTAU;
			++pair_cnts.tau_threshold;
		}


		{
			++pair_cnts.output;
			gene_flags[gi1] |= eOUTPUT;
			gene_flags[gi2] |= eOUTPUT;
		}


		cout	<< name1		<< "\t"
				<< name2		<< "\t";
		switch (output_mode)
		{
		case edS:
			cout << dndsdata.dS << "\n";
			break;
		case edS_COVERAGE:
			cout << dndsdata.dS << "\t"
				 << index_of_fragment << "\n";
			break;
		case edN:
			cout << dndsdata.dN << "\n";
			break;
		default:
			cout << dndsdata << "\t"
				 << aligndata << "\n";
		}
	}
};





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_program_args args;

	try
	{

	debug_func_cerr(5);

	// don't use C streams
	std::ios::sync_with_stdio(false);


	// get command line arguments
	if (!process_command_line_options(argc, argv, args))
		return 1;



	//
	// filter alignment data
	//
	args.VERBOSE() << "\tReading sequence pair data from [" + 
						(args.file_input.length() ? 
						 get_filename(args.file_input) : 
						 "STDIN") + "]\n";
	t_threshold_counts pair_cnts("Alignment pairs", args);
	do_parse_pairs parse_pairs(pair_cnts, args,
								args.comment_letters);
	seq_pairs::read_data(parse_pairs, args.file_input, args.VERBOSE());


	// count genes
	t_threshold_counts gene_cnts("identifiers", args);
	count_pairs_over_threshold(gene_cnts, parse_pairs.gene_flags);
//	filter_kaks_data_from_file(pair_cnts, gene_cnts, args);

	// write summary
	args.VERBOSE() << pair_cnts << "\n" << "\n";
	args.VERBOSE() << gene_cnts << "\n" << "\n";
	
	// write to file
	if (!args.results_summary.empty())
	{
		ofstream ofs(args.results_summary.c_str());
		if (!ofs)
			throw runtime_error("Could not open " + args.results_summary);
		ofs << pair_cnts << "\n" << "\n";
		ofs << gene_cnts << "\n" << "\n";
	}

	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

